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Abstract 

We introduce a new equation describing epitaxial growth processes. This equation is derived 
from a simple variational geometric principle and it has a straightforward interpretation in terms 
of continuum and microscopic physics. It is also able to reproduce the critical behavior already 
observed, mound formation and mass conservation, but however does not fit a divergence form 
as the most commonly spread models of conserved surface growth. This formulation allows us 
to connect the results of the dynamic renormalization group analysis with intuitive geometric 
principles, whose generic character may well allow them to describe surface growth and other 
phenomena in different areas of physics. 
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Variational principles can be found everywhere in physics. Classical mechanics, optics, 
electrodynamics, and general relativity are examples of fundamental theories that are de- 
scribed by one such principle. These principles are not only a mathematically elegant way of 
describing the theory, but they are also supposed to express the fundamental physical con- 
tent of it. Quantum mechanics and quantum field theories, the most precise theories ever 
developed in any scientific discipline, can be expressed as quantizations (e. g. by means 
of path integrals) of the corresponding classical variational principles. The main role that 
variational principles have played in a vast range of physical theories makes them desirable 
when building new theoretical frameworks. 

Non-equilibrium statistical mechanics is an area that has undergone an important growth 
in the last decades. The large number of phenomenologies found in non-equilibrium systems 
makes it difficult to develop unifying principles as in other fields of physics. Among all the 
areas of non-equilibrium statistical mechanics, few of them have the relevance and practical 
importance of surface growth. Together with important technological processes like thin film 
deposition, it describes many phenomena driven by growing interfaces in physics, biology, 
and other scientific fields [1] . The goal of the present work is to develop a variational theory 
of surface growth and, in particular, in the case of epitaxial growth. Our theory will be based 
on a simple geometric principle and it will, at the same time, express plausible microscopic 
processes (which describe the adatoms dynamics) and the observed macroscopic behavior. 

Epitaxial growth is characterized by the deposition of new material on existing layers of 
the same material under high vacuum conditions. The mathematical description uses the 
function h = h{x,y,t), which describes the height of the growing interface in the spatial 
point (x, y) at time t. Although this theoretical framework can be extended to any spatial 
dimension d, we will concentrate here on the physical situation d = 2. A basic assumption 
is the no overhang approximation, which corresponds with the possibility of parameterizing 
the interface as a Monge patch. The macroscopic description of the growing interface is 
given by a stochastic partial differential equation (SPDE) which is usually postulated using 
phenomenological arguments. Examples of such theories are given by the well known SPDEs 
named after Kardar, Parisi, and Zhang (KPZ) [2], Edwards and Wilkinson (EW), MuUins 
and Herring (MH) [l\ and Villain, Lai, and Das Sarma (VLDS) [3]. The last three are of 
the form 

dth = -V-i + F + r]{^,t), (1) 
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where the constant deposition rate F will be absorbed in the definition of the function height 
by means of a Galilean transformation of the reference frame /;,(x, t) — > h{'x,t) + Ft for the 
reminder of this work. The spatial variable x = (x, y) and the noise term 77 is a Gaussian 
distributed random variable that takes into account the thermal fluctuations. It has zero 
mean and its correlation is given by (77(x, t)?7(x', t')) = D5(x — x.')6(t — t'), where D is the 
noise intensity. Eq. ([1]) expresses the conservation of the current j, while the noise is not 
conserved. During years, the most general equation of this type, up to leading order, has 
been supposed to be 

dth = z/2V'/i + A22 V'(V/i)' + XnV{Vhf - v^V^h + r]{^, t), (2) 

the discussion being restricted, mainly, to which terms need to be included, and which do 
not. The aforementioned SPDEs are obtained when we set to zero some of these coefficients, 
say, EW implies A22 = A13 = z/4 = 0, MH corresponds to z/2 = A22 = A13 = 0, and VLDS is 
recovered when i>2 = A13 = 0. However, as we will show in a moment, the divergence form 
Eq. ([1]) and the leading order Eq. ([2]) are not the most general ways of expressing the drift 
term of a conserved growth mechanism. 

Macroscopic equations of surface growth, as Eq. ([2]), have been usually derived according 
to phenomenological arguments. An exception to this has been the recent derivations directly 
from discrete microscopic models, using suitable modifications of the van Kampen system 
size expansion j4|. Herein we will use a different approach: we will consider the variational 
formulation of surface growth following the seminal idea introduced in j^. In order to 
proceed with our derivation, we will assume that the height function obeys a gradient flow 
equation 

dh 6V ^ 

where we have added the noise term, what could be thought as a Parisi-Wu stochastic 
quantization of the corresponding classical variational principle. The functional V denotes a 
potential that is pursued to be minimized during the temporal evolution of h. This potential 
describes the microscopic properties of the interface and of the adatom interactions and, at 
large enough scales, we assume that it can be expressed as a function of the surface mean 
curvature only: 

V = [ fiH)^d^, (4) 
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where H denotes the mean curvature, g the determinant of the surface metric tensor, and 
/ is an unknown function of H. We will further assume that this function can be expanded 
in a power series 

f(H) = K + K^H+^H' + ... , (5) 

of which only the zeroth, first, and second order terms will be of relevance at large scales. 
The result of the minimization of the potential (jlj) leads to the SPDE 

dth = KV^h + K^ [{d,,h)idyyh) - {d,yh)^] - K^V^ h + T] , t) , (6) 

to leading order in the small gradient expansion, which assumes \Vh\ <^ 1 [6]. Note that 
this equation possesses the same linear terms as the generic equation ([2]), but however the 
nonlinearity is rather different; it is a well known nonlinear differential operator referred 
to as the Monge-Ampere operator. The properties of this equation are quite remarkable. 
First of all, the drift is not of divergence form, but it anyway expresses mass conservation. 
This can be easily seen, for instance, assuming that the growth domain is the rectangle 
[—L, L] X [—L, L] with no flux boundary conditions. Then, integrating Eq. ([6]) over the 
whole domain the linear terms vanish immediately and the nonlinear term after integrating 
by parts twice. The conclusion is that the integral of the function h over the whole domain, 
this is the total mass, is conserved over time. 

Apart from the conservation of mass, the scaling properties of this equation are also the 
desired ones. If we introduce the usual scaling x e'x, y — e'y, t eH, and h e'"/i 
with / > 0, we find that scale invariance is found for the EW exponents a = and z = 2, 
and so (3 = a/z = 0, because the terms proportional to Ki and K2 become irrelevant in 
the renormalization group sense. To see what are the intermediate dynamics before the EW 
fixed point dominates we will neglect the term proportional to K and study the equation 

dth = Ki [id,.h)idyyh) - id^^yhf] - K^W^h + r/(x, t), (7) 

which can be thought as the counterpart of the VLDS equation. Applying the standard 
dynamic renormalization group analysis j3] to this SPDE we arrive at the following renor- 
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malization group flow equations at one loop order 
rIK 

^ = K,{a-A + z), (9) 
^ = D{z-2-2a), (10) 

where 6 is the angle formed by the wavevector and the abscise axis, which parameterizes the 
surface in radial coordinates. From here we can straightforwardly read that scale invariance 
is reached for the values of the critical exponents a = 2/3, z = 10/3, and f] = 1/5, this 
is, Eq. ([7]) is in the VLDS universality class. The renormalization group analysis clearly 
shows that the term proportional to K2 is irrelevant (up to one loop order), and that the 
interface dynamics is dominated by the Ki term in the long time limit, a result identical 
to the one obtained for the VLDS equation j^. There is, however, one important difference 
among both equations respecting to the renormalization of K2. For the VLDS equation, the 



renormalization group flow ([8 
that of the ID KPZ equation 



approaches a positive fixed point, a situation reminiscent to 
. In the KPZ situation, however, the interface dynamics is a 
nontrivial combination of diffusion and nonlinearity, while in the VLDS one the nonlinearity 
dominates. In the present case, except for 6 = 7r/2,37r/2, Eq. ([8]) approaches zero with an 
increasing velocity till eventually blowing up dK2/dl —00 when K2 0. So this point is 
the end of the classical evolution of the differential equation and, in our current view, the 
end of the physical evolution too. Anyway, the point K2 = is a branching point, and one 
could consider extensions of the solution to a different branch, negative or complex (allowing 
negative branches will rend the theory unstable and one would have to continue expansion ([5]) 
in order to recover stability \^]). But we are not going to consider such extensions as physical 
within this context, li 6 = 7r/2,37r/2, then K2{1) simply falls exponentially to zero. These 
results explain how the nonlinearity becomes dominant and the diffusion irrelevant in the 
hydrodynamic limit. One may wonder about the 6 dependence in Eq. (jH]). It means that 
the renormalization of diffusion constant K2 depends on the nonlinearity, which is actually 
the small gradient approximation of the surface Gaussian curvature (see below). Diffusion is 
only affected by curvature along the direction of diffusion, and thus the explicit dependence 
on 6 after fixing the coordinate system. Of course, the result becomes independent of this 
angle asymptotically in the scale /, as we have already seen. 
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So we see that if Eq. ([6]) describes an interface that initially fluctuates following MH 
dynamics (this is, a = 1, z = 4, and jS = 1/4), then it will switch to the VLDS universality 
class for intermediate times, and finally, in the long time, will behave as the EW equation. 
This same transient of scaling behavior was found in [sl by means of a formidable renormal- 
ization group calculation of the 2D activated surface diffusion model, except for the final EW 
behavior. These authors also find the MH initial condition for the renormalization group 
flow equations, preceded in certain temperature ranges by a short transient characterized 
by conserved MH dynamics. Of course, this initial transient cannot be found in our theory, 
since we are not considering conserved fluctuations. One obvious improvement would be to 
include both conserved and nonconserved noise in a generalization of the Parisi-Wu quanti- 
zation ([3]). The asymptotic EW dynamics was found in the 2D EW discrete model, in which 
a special type of adatom diffusion is considered j^. Presumably, combining both models a 
similar output to ours could be found. Herein we have seen that such evolutions have a clear 
geometrical meaning: short times are characterized by the minimization of the square of the 
mean curvature, during intermediate times the mean curvature itself is minimized, while in 
the long time limit the zeroth power of the curvature, which corresponds to the surface area, 
is minimized. These dynamical behaviors are in correspondence with the MH, VLDS, and 
EW universality classes respectively. Furthermore, this theoretical framework allows us to 
continue the series expansion in terms of powers of the curvature to include cubic or higher 
orders, and so to describe shorter spatiotemporal scales, at least in principle. Practically it 
might difficult to perform this program up to an arbitrary order in the expansion due to the 
possible presence of nonuniversal features when the spatial scales become of the same order 
of magnitude of the lattice spacing. 

Another characteristic feature that has been observed in both numerical simulations and 
experiments of epitaxial growth is mound formation [4]. Eq. ([6]) is able to reproduce this 
feature too. The terms proportional to K and K2 are stabilizing, and push the surface 
towards planar profiles. On the other hand, the nonlinearity proportional to Ki favors 
mound formation. The Gaussian curvature of the surface is given by 

id^^h){dyyh) - {d^yhf 

= iTWW ' ^ ^ 

which is of the same sign as the Ki term. This shows that those parts of the interface 
provided with a positive Gaussian curvature grow, while those provided with a negative 
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Gaussian curvature decrease in time. This indicates how mounds alone might form, contrary 
to hnear instabihties that favor both negative and positive curvatures, and thus mounds 
and valleys formation. To further clarify how mounds appear, we will consider the evolution 
equation for just the nonlinear term 

dth = K, [{d^^h) (dyyh) - id^yhy] . (12) 

By means of a variables separation h{x, y, t) = a{t)b{x, y) we get 

C, (13) 



a' _ {d^xb){dyyb) - {d^yhf 



K^a? b 
for some constant C. The solution for the temporal variable reads 

a{t) = {[a{0)]-'-CK,ty\ (14) 

showing that for C < it decays in time, while for C > it grows unboundedly till it blows 
up in finite time. Of course, in the real situation blow up will not occur due to the under- 
lying lattice structure. Because mass is conserved, the blow up appears as the consequence 
of the collapse of a microscopic structure. Such a collapse is precluded by the existence of a 
finite lattice spacing, and it is affected by the regularizing effects of the K and K2 terms. By 
Eq. (JT3II we find, as expected, that the sign of the constant C is the same that the sign of the 
Gaussian curvature of the surface. We only need to show that there exist nontrivial positive 
solutions to the equation {dxxb){dyyb) — {dxybY = Cb, for some positive constant C. The 
existence of such a solution was rigorously proven in [9], what opens the possibility of this 
mechanism of mound formation. We do not know if there exist nontrivial solutions to this 
equation for some negative values of the constant C, but this fact is physically unimportant, 
as these solutions, provided they exist, will approach the trivial one monotonically in time. 
One important point here is that the nonlinear Monge- Ampere term is able to promote both 
mound formation and VLDS critical dynamics. Both features were found simultaneously 
in heteroepitaxial growth of InP layers [l^, and so this mechanism constitutes a plausible 
explanation of the observed phenomenology. Note that we can consider Eq. ([7]) as a pos- 
sible continuum description of the 2D activated surface diffusion model (which describes 

n 

homoepitaxial growth) or the heteroepitaxial growth experiments [10] . On the other hand, 
the full equation ([6]) can be considered as the description of the epitaxial system after a 
suitable control protocol, incorporating EW diffusion, has been introduced. The resulting 



dynamics will approximate asymptotically the EW universality class, and thus the desired 
the smooth (or more precisely logarithmically rough) interface. 

The variational formulation allows us to derive the stationary probability functional of 
the interface configuration. In the small gradient approximation the potential functional 
reads 



V[/i(x,t)] 



-f (V/l)' + K^KhyKy - ^{V'hf 



dxdy, (15) 



and so, the probability of a given spatial configuration h{x.) in the long time limit is given 
by 

Vs[h{^)] = Af exp {-2D-^V[h{^)]} , (16) 
where the normalization constant is 

AT-i = J Vh{x) exp {-2D-^V[h{x)]} , (17) 

and the functional integral extends to all the functions h{x) satisfying the prescribed bound- 
ary conditions. Another advantage of Eq. ([6]) is its easy interpretation in terms of microscopic 
physical processes, this is, the underlying adatoms dynamics. The term proportional to K 
is simply diffusion, and correspondingly denotes the random walk character of the adatom 
motion on the surface. The term proportional to K2 is microscopically related to fine prop- 
erties of the random walk. Such a term may arise, for instance, when studying corrections 
to the central limit theorem due to a small scale underlying structure. It is clear now why 
it is the first in becoming irrelevant when subject to renormalization group flow. The non- 
linear term, proportional to Ki, expresses optimal mass rearrangement: the new deposited 
adatoms, which form a disordered structure, move towards the new, more favorable fiat 
disposition that minimizes the interface chemical potential. This rearrangement is such that 
the square of the distance among the initial and final adatoms position is minimized 11|. 
The reason why the squared distance, and not the distance itself, is minimized has to do 
with the random dispersal properties of the adatoms. As they move randomly, they propa- 
gate along areas rather than distances. Note that this diffusive process could be identified 
with the driving force generating mounds in the experimental situation described in flo| . 
It is remarkable that a similar process of optimal mass rearrangement has been found in 
atmospheric physics, particularly in semigeostrophic flow |l2|. 

The nonlinear term in Eq. ([6]) was derived from a specific variational principle, but it 
appears in much more general situations. It can actually be derived from an infinite number 
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of inequivalent potential functionals, a fact specially relevant in the context of high energy 



physics 



13|. Apart from its interesting geometrical properties (it is related to both mean and 



Gaussian curvatures of the surface, as we have seen), it also appears in the theory of optimal 



transport 



11| and in this context is related to atmospheric air flow [13] ■ Another interesting 



property of Eq. (jH]) is its connection with the ID KPZ equation: assuming the solution 
form h{x,y,t) = —yu{x,t), we see that it obeys the ID KPZ equation after including y 
dependence in the coefficient of the nonlinear term and setting K2 = 0. Despite the purely 
mathematical nature of this symmetry, the generic character of both equations may reveal 
its physical significance in a field different from surface growth. 

In conclusion, we have derived a stochastic growth equation from a simple variational 
principle. This equation describes the physics of the growing interface at both the macro 
and micro levels, and reproduces mound formation and the known critical exponents. Such 
a complete and unified approach has not been achieved previously, to our knowledge, in the 
epitaxial growth literature. This is so possibly due to the omnipresent assumption of the 
divergence form Eq. ([T]) for models describing epitaxial growth. In fact, this assumption 
has been adopted for years as the only way of assuring the conservation of mass; however, 
as we have shown, mass is conserved in models not belonging to this class. Furthermore, 
Eq. ([6]) can be derived from simple and general principles, what makes it a very generic 
construction of potential influence in a wide range of different areas within physics. As we 
have already mentioned, high energy, atmospheric, condensed matter, and statistical physics 
can be connected to it. Indeed, we are persuaded that the profoundly generic nature of this 
construction will make it appear in many other different SbTQdbS clS well. 
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